home *** CD-ROM | disk | FTP | other *** search
/ Over 1,000 Windows 95 Programs / Over 1000 Windows 95 Programs (Microforum) (Disc 1).iso / 1249 / sunup.t < prev    next >
Encoding:
Text File  |  1997-04-18  |  9.3 KB  |  442 lines
  1. %
  2. % "sunup.t" calculates the times of sunrise and sunset
  3. % on any date, accurate to the minute within several
  4. % centuries of the present.  It correctly describes what
  5. % happens in the arctic and antarctic regions, where
  6. % the Sun may not rise or set on a given date.  
  7. % 
  8. % This program was adapted from a BASIC program in
  9. % Sky & Telescope magazine, August 1994, page 84,
  10. % written by Roger W. Sinnott.
  11. %
  12. % Enter north latitudes positive, west longitudes negative.  
  13. %
  14. % For the time zone, enter the number of hours west of Greenwich 
  15. % (e.g., 5 for EST, 4 for EDT).  
  16. %
  17. %   Sample program for the T Interpreter by:
  18. %
  19. %   Stephen R. Schmitt
  20. %   962 Depot Road
  21. %   Boxborough, MA 01719
  22. %
  23.  
  24. const PI : real := 2 * arcsin( 1 )
  25. const DR : real := PI / 180
  26. const K1 : real := 15 * DR * 1.0027379
  27.  
  28. var Sunrise, Sunset : boolean := false
  29.  
  30. type rvector : array[3] of real
  31.  
  32. program
  33.  
  34.     var lat, lon, tz, ph : real
  35.     var jd, ct, t0, ra0, ra1, dec0, dec1 : real
  36.     var k, zone : int
  37.     var ra, dec, v : rvector
  38.  
  39.     prompt "Latitude, +N, -S (deg):"
  40.     get lat
  41.     prompt "Longitude, +E, -W (deg):"
  42.     get lon
  43.     prompt "Time zone +West of GMT, -East (hrs):"
  44.     get zone
  45.  
  46.     jd := julian_day - 2451545  % Julian day relative to Jan 1.5, 2000
  47.     show_input( lat, lon, zone )
  48.     
  49.     lon := lon / 360
  50.     tz := zone / 24
  51.     ct := jd / 36525 + 1        % centuries since 1900.0
  52.     
  53.     t0 := lst( lon, jd, tz )
  54.     
  55.     % get sun's position
  56.     jd := jd + tz
  57.     sun( ra0, dec0, jd, ct )
  58.     jd := jd + 1
  59.     sun( ra1, dec1, jd, ct )
  60.  
  61.     if ra1 < ra0 then 
  62.         ra1 := ra1 + 2 * PI
  63.     end if
  64.  
  65.     ra[0]  := ra0
  66.     dec[0] := dec0
  67.     
  68.     for k := 0...23 do
  69.  
  70.         ph := ( k + 1 ) / 24
  71.         
  72.         ra[2]  := ra0  + ph * ( ra1  - ra0 )
  73.         dec[2] := dec0 + ph * ( dec1 - dec0 )
  74.     
  75.         test( k, zone, t0, lat, ra, dec, v )
  76.         
  77.         ra[0]  := ra[2]
  78.         dec[0] := dec[2]
  79.         
  80.     end for
  81.     
  82.     %  special message?
  83.     if ( not Sunrise ) and ( not Sunset ) then
  84.  
  85.         if v[2] < 0 then
  86.             put "Sun down all day"
  87.         else
  88.             put "Sun up all day"
  89.         end if
  90.  
  91.     else
  92.     
  93.         if not Sunrise then
  94.             put "No sunrise this date"
  95.         elsif not Sunset then
  96.             put "No sunset this date"
  97.         end if
  98.  
  99.     end if
  100.  
  101. end program
  102.  
  103. %
  104. % display latitude, longitude and time zone
  105. %
  106. procedure show_input( lat, lon : real, zone : int )
  107.  
  108.     var deg, min : int
  109.     
  110.     if sgn( zone ) = sgn( lon ) and zone ~= 0 then
  111.         put "WARNING: time zone and longitude are incompatible"
  112.     end if
  113.  
  114.     if abs( zone ) > 12 then
  115.         put "WARNING: invalid time zone"
  116.     end if        
  117.  
  118.     if abs( lat ) > 90 then
  119.         put "WARNING: invalid latitude"
  120.     end if
  121.  
  122.     if abs( lon ) > 180 then
  123.         put "WARNING: invalid longitude"
  124.     end if
  125.  
  126.     if lat > 0.0 then
  127.         deg := floor( lat )
  128.         min := floor( 60 * ( lat - deg ) )
  129.         put deg:4, ":",min:2, " N "...
  130.     else
  131.         deg := floor( -lat )
  132.         min := floor( 60 * ( -lat - deg ) )
  133.         put deg:4, ":",min:2, " S "...
  134.     end if
  135.  
  136.     if lon > 0.0 then
  137.         deg := floor( lon )
  138.         min := floor( 60 * ( lon - deg ) )
  139.         put deg:4, ":",min:2, " E"
  140.     else
  141.         deg := floor( -lon )
  142.         min := floor( 60 * ( -lon - deg ) )
  143.         put deg:4, ":",min:2, " W"
  144.     end if
  145.     
  146. end procedure
  147.  
  148. %
  149. % LST at 0h zone time
  150. %
  151. function lst( lon, jd, z : real ) : real
  152.  
  153.     var s, t0 : real
  154.     
  155.     s := 24110.5 + 8640184.812999999 * jd / 36525
  156.     s := s + 86636.6 * z + 86400 * lon
  157.     
  158.     s := s / 86400
  159.     s := s - floor( s )
  160.  
  161.     t0 := s * 360 * DR
  162.  
  163.     return t0
  164.     
  165. end function
  166.  
  167. % 
  168. % test an hour for an event
  169. %
  170. procedure test( k, zone : int, 
  171.                 t0, lat : real,
  172.                 var ra, dec, v : rvector )
  173.  
  174.     var ha : rvector
  175.     var a, b, c, d, e, s, z : real
  176.     var hr, min, time : real
  177.     var az, dz, hz, nz : real
  178.     var zchar : string
  179.     var zlet : string := "MLKIHGFEDCBAZNOPQRSTUVWXY"
  180.     label test_exit :
  181.     
  182.     ha[0] := t0 - ra[0] + k * K1 
  183.     ha[2] := t0 - ra[2] + k * K1 + K1 
  184.  
  185.     ha[1]  := ( ha[2]  + ha[0]  ) / 2       % hour angle at half hour
  186.     dec[1] := ( dec[2] + dec[0] ) / 2       % declination at half hour
  187.     
  188.     s := sin( lat * DR )
  189.     c := cos( lat * DR )
  190.     z := cos( 90.833 * DR )
  191.  
  192.     if k <= 0 then
  193.         v[0] := s * sin( dec[0] ) + c * cos( dec[0] ) * cos( ha[0] ) - z
  194.     end if
  195.  
  196.     v[2] := s * sin( dec[2] ) + c * cos( dec[2] ) * cos( ha[2] ) - z
  197.     
  198.     if sgn( v[0] ) = sgn( v[2] ) then
  199.         goto test_exit
  200.     end if
  201.     
  202.     v[1] := s * sin( dec[1] ) + c * cos( dec[1] ) * cos( ha[1] ) - z
  203.     
  204.     a :=  2 * v[0] - 4 * v[1] + 2 * v[2] 
  205.     b := -3 * v[0] + 4 * v[1] - v[2]
  206.     
  207.     d := b * b - 4 * a * v[0]
  208.  
  209.     if d < 0 then
  210.         goto test_exit
  211.     end if
  212.     
  213.     d := sqrt( d )
  214.     
  215.     if v[0] < 0 and v[2] > 0 then
  216.         put "Sunrise at: "...
  217.         Sunrise := true
  218.     end if
  219.     
  220.     if v[0] > 0 and v[2] < 0 then
  221.         put "Sunset at:  "...
  222.         Sunset := true
  223.     end if
  224.     
  225.     e := ( -b + d ) / ( 2 * a )
  226.     
  227.     if e > 1 or e < 0 then
  228.         e := ( -b - d ) / ( 2 * a )
  229.     end if
  230.     
  231.     % time of an event
  232.  
  233.     time := k + e + 1 / 120
  234.     
  235.     hr := floor( time )
  236.     min := floor( ( time - hr ) * 60 )
  237.  
  238.     zchar := "    "
  239.     if abs( zone ) <= 12 then
  240.         zchar[1] := zlet[zone+12]
  241.     end if
  242.  
  243.     put hr:2, ":", min:2, zchar...
  244.  
  245.     % azimuth of the sun at the event
  246.  
  247.     hz := ha[0] + e * ( ha[2] - ha[0] )
  248.     nz := -cos( dec[1] ) * sin( hz )
  249.     dz := c * sin( dec[1] ) - s * cos( dec[1] ) * cos( hz )
  250.     
  251.     az := arctan( nz / dz ) / DR
  252.     
  253.     if dz < 0 then
  254.         az := az + 180
  255.     end if
  256.     
  257.     if az < 0 then
  258.         az := az + 360
  259.     end if
  260.     
  261.     if az > 360 then
  262.         az := az - 360
  263.     end if
  264.     
  265.     put "azimuth: ", az:5:1
  266.  
  267. test_exit:
  268.  
  269.     v[0] := v[2]
  270.  
  271. end procedure
  272.  
  273. % 
  274. % sun's position using fundamental arguments 
  275. % (Van Flandern & Pulkkinen, 1979)
  276. %
  277. procedure sun( var ra, dec : real, jd, ct : real )
  278.  
  279.     var g, lo, s, u, v, w : real
  280.     
  281.     lo := 0.779072 + 0.00273790931 * jd
  282.     lo := lo - floor( lo )
  283.     lo := lo * 2 * PI
  284.  
  285.     g := 0.993126 + 0.0027377785 * jd
  286.     g := g - floor( g )
  287.     g := g * 2 * PI
  288.     
  289.     v := 0.39785 * sin( lo )
  290.     v := v - 0.01 * sin( lo - g )
  291.     v := v + 0.00333 * sin( lo + g )
  292.     v := v - 0.00021 * ct * sin( lo )
  293.     
  294.     u := 1 - 0.03349 * cos( g )
  295.     u := u - 0.00014 * cos( 2 * lo )
  296.     u := u + 0.00008 * cos( lo )
  297.     
  298.     w := -0.0001 - 0.04129 * sin( 2 * lo )
  299.     w := w + 0.03211 * sin( g )
  300.     w := w + 0.00104 * sin( 2 * lo - g )
  301.     w := w - 0.00035 * sin( 2 * lo + g )
  302.     w := w - 0.00008 * ct * sin( g )
  303.     
  304.     % compute sun's right ascension and declination
  305.     
  306.     s := w / sqrt( u - v * v )
  307.     ra := lo + arctan( s / sqrt( 1 - s * s ) )
  308.     s := v / sqrt( u )
  309.     dec := arctan( s / sqrt( 1 - s * s ) )
  310.  
  311. end procedure
  312.  
  313. %
  314. % determine Julian day from calendar date
  315. % ref: Jean Meeus, "Astronomical Algorithms", Willmann-Bell, 1991
  316. %
  317. function julian_day : real
  318.  
  319.     var a, b, day, jd : real
  320.     var month, year : int
  321.     var gregorian : boolean
  322.  
  323.     prompt "Year:"
  324.     get year
  325.     prompt "Month:"
  326.     get month
  327.     prompt "Day:"
  328.     get day
  329.  
  330.     date( floor( day ), month, year )
  331.     
  332.     if year < 1583 then
  333.         gregorian := false
  334.     else
  335.         gregorian := true
  336.     end if
  337.     
  338.     if month = 1 or month = 2 then
  339.         year := year - 1
  340.         month := month + 12
  341.     end if
  342.  
  343.     a := floor( year / 100 )
  344.     if gregorian then
  345.         b := 2 - a + floor( a / 4 )
  346.     else
  347.         b := 0.0
  348.     end if
  349.  
  350.     jd := floor( 365.25 * ( year + 4716 ) ) + 
  351.           floor( 30.6001 * ( month + 1 ) ) + 
  352.           day + b - 1524.5
  353.     
  354.     return jd
  355.     
  356. end function
  357.  
  358. %
  359. % display the date
  360. %
  361. procedure date( d, m, y : int )
  362.  
  363.     var day, month : array[12] of string
  364.     var mi, yi, i : int
  365.  
  366.     month[0]  := "January"
  367.     month[1]  := "February"
  368.     month[2]  := "March"
  369.     month[3]  := "April"
  370.     month[4]  := "May"
  371.     month[5]  := "June"
  372.     month[6]  := "July"
  373.     month[7]  := "August"
  374.     month[8]  := "September"
  375.     month[9]  := "October"
  376.     month[10] := "November"
  377.     month[11] := "December"
  378.  
  379.     day[0] := "Monday"
  380.     day[1] := "Tuesday"
  381.     day[2] := "Wednesday"
  382.     day[3] := "Thursday"
  383.     day[4] := "Friday"
  384.     day[5] := "Saturday"
  385.     day[6] := "Sunday"
  386.  
  387.     mi := m
  388.     yi := y
  389.  
  390.     if y < 1752 then
  391.  
  392.         put month[mi-1], " ", d, ", ", yi
  393.  
  394.     else
  395.  
  396.         if m = 1 or m = 2 then
  397.  
  398.             m := m + 12
  399.             y := y - 1
  400.  
  401.         end if
  402.  
  403.         i := d + 2*m + 3*(m+1) div 5 + y + y div 4 - y div 100 + y div 400
  404.         i := i mod 7
  405.  
  406.         put day[i], ", ", month[mi-1], " ", d, ", ", yi
  407.  
  408.     end if
  409.  
  410. end procedure
  411.  
  412. %
  413. % returns value for sign of argument
  414. %
  415. function sgn( x : real ) : int
  416.  
  417.     var rv : int
  418.     
  419.     if x > 0.0 then
  420.         rv :=  1
  421.     elsif x < 0.0 then
  422.         rv := -1
  423.     else
  424.         rv :=  0
  425.     end if
  426.  
  427.     return rv
  428.  
  429. end function
  430.  
  431. %
  432. % absolute value of argument
  433. %
  434. function abs( x : real ) : real
  435.  
  436.     if x < 0.0 then
  437.         x := -x
  438.     end if
  439.  
  440.     return x
  441.  
  442. end function